Subdiffusion of nonlinear waves in quasiperiodic potentials 
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We study the time evolution of wave packets in one-dimensional quasiperiodic lattices which 
localize linear waves. Nonlinearity (related to two-body interactions) has a destructive effect on 
localization, as recently observed for interacting atomic condensates [Phys. Rev. Lett. 106, 230403 
(2011)]. We extend the analysis of the characteristics of the subdiffusive dynamics to large temporal 
and spatial scales. Our results for the second moment m,2 consistently reveal an asymptotic 7712 ~ 
t 1 ^ 3 and intermediate ni2 ~ t 1 ^ 2 laws. At variance to purely random systems [Europhys. Lett. 91, 
30001 (2010)], the fractal gap structure of the linear wave spectrum strongly favors intermediate 
self-trapping events. Our findings give a new dimension to the theory of wave packet spreading in 
localizing environments. 
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I. INTRODUCTION 

In one dimension, a wave packet of noninteracting par- 
ticles subject to a random potential does not diffuse be- 
cause of Anderson localization, due to the exponential 
localization of the eigenstates of the underlying Hamil- 
tonian [J 0]. Instead the presence of interactions is 
expected to act against localization, though the actual 
mechanism may be highly nontrivial and may depend on 
the type of disorder and interaction. This is a problem 
of fundamental importance for many systems in different 
contexts. From the theoretical side, it has been largely 
studied by using discrete lattices. The interaction is of- 
ten included by means of mean- field theories, where it 
enters in the form of a nonlinear local term in a nonlin- 
ear Schrodinger-like equation (3rll3j. 

Numerical simulations of wave packets propagating in 
a random potential - with the interaction included via 
a nonlinear mean-field term - showed that the presence 
of interaction indeed destroys localization and leads to 
a subdiffusive growth of the second moment of the wave 
packet in time as t 1 Eol4l3 |. In particular it was 
predicted that at large t, the coefficient 7 should con- 
verge to 1/3 in a regime of so-called "weak chaos", as 
opposed to normal diffusion where 7 = 1 and the wave 
packet width grows as \ft. A transient regime of "strong 
chaos" was also identified, where 7 = 1/2 [3, El El- The 
occurrence of these different regimes can be predicted by 
comparing the nonlinear frequency shift introduced by 
the expanding wave packet to the typical energy scales 
of linear spectra of random models. 

Exponential localization for noninteracting quantum 
particles (or linear waves) can be found also in systems 
which are not truly disordered. An example is provided 
by quasiperiodic potentials, which are of great interest by 
themselves 0, EH ■ These systems can be considered to 



lie between the two extreme cases of a perfectly periodic 
system and a pure random potential. By tuning the pa- 
rameters of a quasiperiodic system, the localization prop- 
erties can change dramatically - from having all extended 
states to all localized states. In recent years, exponential 
localization has been observed with light propagating in 
quasiperiodic photonic lattices [Hij , as well as with ultra- 
cold atoms propagating in a bichromatic optical lattice 
fl7j . Notably in both cases the inclusion of interaction is 
experimentally feasible, by using a Kerr medium for light 
and tuning the scattering length by means of a suitable 
magnetic field for atoms. 

Numerical simulations studying nonlinear dynamics of 
wave packets have been performed also in the case of 
quasiperiodic systems |18l - [21| . In particular for expo- 
nentially localized linear waves, nonlinearity yields subd- 
iffusive spreading of wave packets as well [20(. However, 
there are clear indications that the coefficient 7, at least 
at finite spreading times, is significantly larger than the 
one observed in random systems. Nonlinear effects have 
been also studied in experiments using ultracold atoms 
and light propagating in photonic lattices. In both cases 
it has been shown that nonlinearity acts against localiza- 
tion mm. 

The purpose of the present work is to clarify the de- 
tails of the spreading mechanism leading to the destruc- 
tion of the localization in quasiperiodic systems, and to 
address differences and similarities between quasiperiodic 
and purely random potentials. We extend and refine pre- 
vious numerical investigations by pushing the simulations 
to much longer times, thus allowing for the identification 
of the strong and weak chaos regimes in quasiperiodic sys- 
tems and compare the situation with known properties of 
purely random systems. For this purpose, we use two dif- 
ferent models; namely, a discrete nonlinear Schrddinger 
equation (DNLS) and a quasiperiodic version of the quar- 
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tic Klein-Gordon (KG) lattice model. 

A key result of the present work is that a regime of 
weak chaos is indeed observed in the long time spread- 
ing of nonlinear wave packets propagating in quasiperi- 
odic systems; in particular we find that the asymptotic 
value of the spreading coefficient 7 is 1/3 as in purely 
random systems, thus showing that this behaviour is 
rather general and model independent. Another similar- 
ity with purely random systems is the occurrence of self- 
trapping: when the nonlinear interaction is large enough 
to shift the mode frequencies so strongly that they are 
tuned out of resonance with all nonexcited neighbour- 
ing modes, apart of the wave packet remains spatially 
localized [H U Hq| • However as opposed to the random 
system, in the quasiperiodic case partial self-trapping is 
also possible for weaker nonlinearities. This is due to the 
complexity of the linear wave spectrum which exhibits 
a fractal gap structure of sub-bands. Self-trapping gives 
rise to transient spreading regimes characterized by an 
intermediate large exponent 7; we call this effect "over- 
shooting". Finally, we have also observed signatures of 
strong chaos, but detection of this regime is difficult in 
quasiperiodic systems, since it is often masked by over- 
shooting and partial self-trapping, which occur on the 
same temporal scales. 



II. MODELS 

We consider two different models: the first is the 
one-dimensional (ID) quasiperiodic discrete nonlinear 
Schrodinger equation (DNLS), defined by the Hamilto- 
nian 



DNLS — 
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(1) 

where j labels the lattice sites and V 3 ■ = A cos(27raj + if). 
The quantity IV^'I 2 gives the probability to find a particle 
at site j. The first term in Eq. (fTJ describes the hop- 
ping between nearest neighbouring sites, the second term 
describes the quasiperiodic on-site energy, while the last 
term represents the mean-field interaction energy and in- 
troduces the nonlinearity. The key parameters of this 
Hamiltonian are the strength of the quasiperiodic poten- 
tial A (for A = the lattice is periodic with period 1), 
the strength of the nonlinearity (3 (for (3 — the particles 
are non- interacting) , and the irrational number a which 
causes the underlying potential to be quasiperiodic. In 
fact, when a is irrational the cosine adds a second pe- 
riodicity which is incommensurate with respect to the 
underlying periodicity given by the discreteness of the 
system. Let us note that, without any loss of generality, 
one can always choose —0.5 < a < 0.5 since Eq. (TT|) and 
the subsequent Eq. ([3]) are invariant under a shift of a by 
an integer number. Choosing the value of a in this in- 
terval has the additional advantage that the wavenumber 



of the underlying discrete lattice. As a convenient choice 
for this work, we use a — (VE— 3)/2 The phase <p is 
a phase shift between the two lattices. The equations of 
motion associated with Eq. (fT]) are idipj/dt = dH/dipj, 
or 



ftk. 
' dt 



(2) 



The above set of equations conserve the energy of Eq. ((T|) , 
as well as the total norm S = Y]j |^j| 2 of the initial wave 
packet (we always assume S = 1). The set can be used for 
a wide range of applications, including ultracold atoms 
expanding in optical lattices Ut I [2 C l | 24l [25l | and light 



propagating in photonic lattices [l6ll26j. 

The second model is a quasiperiodic version of the 
quartic Klein-Gordon (KG) lattice, given by 



H 



KG 



2 ^ 
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where uj and pj are the generalized coordinates and mo- 
menta on the site j and V 3 ■ = 1 + (1/2) cos(27ror/ + (p). 
The energy associated with lattice site j is 

_p 2 V 3 u 2 u] (u J+1 - Uj ) 2 (uj^-Uj) 2 
tj ~ 2 + 2 + 4 + 8A + 8A ' 1 > 

The equations of motion are generated by d 2 Uj/dt 2 = 
—dH/diij, yielding 



dt 2 



3 = -Vjrij 



— (u j+ i +%•_!- 2uj). (5) 



associated to k 



2na, rests in the Brillouin zone 



This set of equations conserve the total energy T~L = 
Ylj £j — only. The KG model has also been exten- 
sively studied, since it can give a simple description of 
the nondissipative dynamics of anharmonic optical lat- 
tice vibrations in molecular crystals [27| . 

The total energy of the system H serves as a control 
parameter of nonlinearity, analogous to (3 for the DNLS 
model. In fact, for small amplitudes the equation of the 
KG chain can be approximately mapped onto a DNLS 
model 28] using (3S « 6XH. Further analytics will be 
discussed only in terms of the DNLS chain, since it is 
then straightforward to project to the KG model. 

For the DNLS model we measure the spreading of the 
wave packet by tracking the quantity rij = | 2 /S, here- 
after named norm density consistently with the notation 
of Refs.fl, 0, H, EH- The key quantities that we use 
to describe the time evolution of the expanding wave 
packet are the second moment m 2 = J2jnj(X — j) 2 , 
(X = n jj) which quantifies the spatial extent of the 
wave packet, and the participation number P = 1/ n 2 
which measures the number of significantly populated 
sites. A combination of these two quantities £ = P 2 /m,2, 
called compactness index, gives a measure of the sparsity 
of a wave packet 7] . For the KG we do exactly the same, 
but replacing the norm density rij with its counterpart 
£j/"H, which is the normalized energy density. 
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III. RELEVANT ENERGY SCALES 

Neglecting the nonlinear term in Eq. ([2]) reduces to an 
eigenvalue problem 

- A v j + i - A ud _i + A coa(2naj + <p)A Vi j = E v A v j . (6) 

Here the index v labels the different normal modes A u j 
and eigenvalues E v . The coefficient 1/(2 A) in Eq. (JS|) was 
chosen so that the linear parts of i?DNLS an d -Hkg would 
correspond to the same eigenvalue problem: the linear 
KG model can then likewise be identically reduced to 
Eq. ©, under the substitution E v = 2A (w^ — 1/A - l), 
where w„ are the corresponding eigenfrequencies. 

Eq. ([6]) is also known as the Aubry- Andre model [Hj]. 
The localization properties of this model are well known 
and extensively studied both analytically and numeri- 
cally 0, QJ, [2(J [H, [H-Hll. A transition occurs from 
an extended regime to a localized regime at A = 2. For 
A < 2 all normal modes A v j are extended over the entire 
lattice, at A = 2 they are critical, while for A > 2 they 
are exponentially localized in the form A u j ~ g-M—iAl 1 ^ 
where j v is the localization center and I — l/ln(A/2) 
is the localization length (notice that it is the same for 
all the modes) [HI]. Since we are interested in the inter- 
play between localization and nonlinearity, we will focus 
exclusively on the regime A > 2. 

In order to quantify the spatial extent of a given eigen- 
state, it is convenient to define a localization volume 

V v = 1 + ^/l2m^, where = ^(X, - j) 2 \A u>j \ 2 

is the second moment of \A„j\ 2 and X v = Y^j j\Av,j\ 2 is 
its center of norm [32j ■ The localization volume is used 
to estimate the number of modes which interact with a 
given mode v. We show its meaning schematically in 
Fig. [TJz. The modes that interact with a given reference 
mode v are those whose center of norm lies in an area 
V v around it. The average localization volume V is then 
found by numerically diagonalizing the linear system, cal- 
culating V v for each eigenmode, and then averaging over 
all eigenmodes. A plot of this quantity as a function of 
the potential strength A is shown in Fig. QJ. 

The spectrum for A > 2 is purely dense-point, charac- 
terized by the presence of an infinite number of gaps and 
bands. A plot of the Aubry-Andre model's spectrum as 
a function of A is shown in Fig. [TJ;. In this figure, one 
clearly sees the presence of two major gaps dividing the 
spectrum in three parts, each of them divided in turn in 
three smaller parts, and so on (33|. We call these portions 
of spectrum separated by the largest gaps "mini-bands" . 
For our purposes, it is enough to consider a division of 
the spectrum in M = 3 or at most in M = 9 mini-bands. 
Smaller mini-bands have vanishingly small effects on the 
time evolution of wave packets. 

Let us introduce two energy scales associated with 
the linear system [f| The first one, A, is the full 

width of the spectrum, defined as the difference be- 
tween the largest and the smallest eigenvalues: A = 
max{E u } — miii{E u }. The second one, d, is the mean 



spacing of eigenvalues within a single mini-band and 
within the range of a localization volume. Let us ex- 
plain how we calculate this quantity. We consider a given 
mini-band and all the eigenstates that lie in it. For each 
eigenstate is, we calculate its localization volume V v and 
then we form the subset of the other eigenstates, {/j,}, 
belonging to the same mini-band and interacting with it, 
namely, those fulfilling the condition \X V — X^\ < V v /2. 
The average number of states in the subset can be esti- 
mated as V/M. Then we calculate the energy spacings 
within this subset. This procedure is repeated for each 
eigenstate in the band and the average gives the mean 
spacing d. 

The number of mini-bands M to be used in the calcu- 
lations of d depends on A. For a given A we choose M in 
such a way that the localization volume V satisfies the 
condition V/M > 2. This implies that, on average, there 
are at least two eigenstates within the subset {/x} that we 
can use to calculate the average energy spacings. We al- 
ways consider A > 2.1; therefore it is enough to divide the 
spectrum at most in nine mini-bands. As A is increased 
the average localization volume of the eigenstates V de- 
creases - therefore at some point we have to consider the 
spectral separation into smaller mini-bands. In practice 
we consider M = 9 mini-bands for 2.1 < A < 2.2, M = 3 
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FIG. 1: (color online) a) Pictorial interpretation of localiza- 
tion volume. A given eigenstate v (black line in the center 
of the box) is assumed to interact only with those eigenstates 
(blue lines) that lie in a region of size V v around his mean 
position. The red lines represent the corresponding on-site 
energies, b) Average localization volume of eigenstates V as 
a function of the potential strength A. c) Eigenenergies E v of 
the linear system obtained from numerical diagonalization of 
Eq. (O, as a function of A. 
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mini-bands for 2.2 < A < 2.75 and just one band (i.e., 
the full spectrum) for A > 2.75. A plot of the energy 
scales A and d as a function of A is shown in Fig. [2l The 
dashed vertical lines represent the values of A where the 
number of mini-bands changes in the calculation of d. 

Similarly to the case of disordered systems 0, 0, [3> 
the scales A and d of the linear spectrum (which are fre- 
quencies in the present setting of nonlinear wave equa- 
tions) must be compared to the frequency shift caused 
by the nonlinearity. Indeed a single oscillator which sat- 
isfies the equation of motion itp = Vip + (3\ip\ 2 ip experi- 
ences a nonlinear frequency shift 8 = fi\ip\ 2 away from 
its linear frequency V. For many oscillators, we can con- 
veniently use the eigenstates of the linear Aubry- Andre 
model as a decomposition basis of the wave function ifij : 
tpj = § v A Vi j. Equation |(3J) can then be rewritten for 
the evolution of the normal mode amplitudes: 

i-^-=E v (p u +ft ^ I v,v 1 ,va,ua < f > ti<Pv a <Pv 3 (7) 

where I v ^ Vl ,v2,v 3 is an overlap integral involving four nor- 
mal modes: 

Iu,vi,u2,u 3 — y ] AvjAu 1 ,jA l/2t jAy 3 j. (8) 

3 

As discussed in Ref.[ll|, one can introduce a norm den- 
sity also in the normal mode space, n u = |</>„| 2 ; as the 
packet spreads and after averaging over many realization, 
this quantity becomes almost identical to the norm den- 
sity rij = \ipj\ 2 and the frequency shift can be expressed 
as 8 ~ fin, where n is a characteristic norm density. In 
the KG model 8 is proportional to the energy density £ 
and, within our formalism, can be obtained by the small 
amplitude mapping. 

When 8 < d, the mode frequencies in a wave packet 
are only weakly shifted, and a small fraction of these 
modes will resonantly and strongly interact with each 
other. Following the terminology of Refs. @, 0, H|, we 
say that this is a regime of weak chaos. Conversely, 
when d < 8 < A, the mode frequencies in a packet are 
strongly shifted and almost all of them will resonantly 
and strongly interact with each other. This is labeled a 
regime of strong chaos. When finally 8 > A, the mode 
frequencies are shifted so strongly that they are tuned 
out of resonance with all nonexcited neighbouring modes. 
An excited mode in this condition may stay localized, i.e., 
self-trapped, for long or even infinite times. The mean- 
ing of these regimes will be further clarified in the next 
section. 



IV. EXPECTED SPREADING REGIMES 

As one can see in Eq. ([7J , the presence of nonlinearity 
in the DNLS model introduces a coupling between eigen- 
states of the underlying linear spectrum. It has already 
been observed numerically and experimentally that this 
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FIG. 2: (color online) Energy scales A (top blue line) and 
d (bottom red line) plotted as a function of the potential 
strength A. The empty (downward) and full (upward) trian- 
gles correspond to the values of S that we have used for the 
simulations with the DNLS model and with the KG model 
respectively. Comparing the nonlinear frequency shift 8 with 
the energy scales A and d one can predict the different spread- 
ing regimes of weak chaos (8 < d), strong chaos (d < 8 < A) 
and self-trapping (8 > A). The separation between the three 
regimes should not be interpreted as a sharp boundary, but 
as a smooth crossover. 



leads to a subdiffusive spreading of wave packets, i.e. its 
second moment grows asymptotically as ni2 ~ t 1 with 
7 < 1 [13, HH- However, a systematic investigation of 
the behaviour of the exponent 7 in different regimes of 
strong and weak chaos, and self-trapping, have not been 
done so far. In this section, we approach this issue by first 
comparing the nonlinear frequency shift 8 = fin with the 
energy scales A and d, in such a way as to introduce the 
different spreading regimes expected to be observed in 
the subsequent numerical simulations. 

Let us consider an initial wave packet with norm den- 
sity n and localization volume L larger than the aver- 
age localization volume of the eigenstates of the linear 
spectrum, L > V. If 8 > A, nonlinearity is so strong 
that all the participating normal modes within the wave 
packet are shifted out of resonance with respect to the 
non-excited neighbourhood; therefore spreading is largely 
suppressed and a significant part of the wave packet re- 
mains self-trapped 34]. If instead 8 < A, we are no 
longer in the self-trapping regime and can distinguish two 
sub-cases: on one hand, when 8 > d, strong chaos is re- 
alized; that is, all the modes in the packet are resonantly 
interacting with each other, thus producing an efficient 
spreading. On the other hand, when 8 < d, weak chaos is 
obtained: only a fraction of modes interact resonantly - 
the localization is still destroyed, but spreading is slower. 

If L < V the estimate of the self-trapping transition 
is done as before, that is by comparing 8 = fin with 
the spectrum width A. If self-trapping is avoided, how- 
ever, the wave packet initially spreads also in absence 
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of nonlinearity, eventually filling the localization volume 
V. Consequently the initial norm density n is reduced 
tons nL/V, due to linear time evolution - the relevant 
nonlinear frequency shift must now be calculated by us- 
ing this reduced density fi. Apart from this detail, which 
originates from the initial dynamics at short times, the 
asymptotic spreading regimes are the same as before. 

Studies performed on random systems have shown that 
the basic mechanism that destroys localization is the 
presence of resonances in mode-mode interaction 
This leads to chaotic dynamics within a part of the wave 
packet, and to a subsequent subdiffusive spreading. Here 
we apply the same theory to the case of quasiperiodic 
systems. For the reader interested in the formal details 
of the theory, we refer to the previous articles @, ls| • 

Let us estimate the number of resonant modes in the 
packet, which is a key quantity that determines the type 
of spreading behaviour. According to Eq. (O, due to 
nonlinearity, the evolution of a given normal mode is af- 
fected by any three (triplet) modes. The coupling is the 
largest if the triplet modes have large amplitudes and if 
the overlap integrals are large, i.e., if the triplet modes are 
close enough in space to the given normal mode. Some of 
these triplet modes may affect the dynamics of the cho- 
sen mode v strongly, some weakly. To distinguish these 
triplet groups, we follow |j[ and apply perturbation the- 
ory to first order in fi. It follows that the amplitude of 
a normal mode v inside the wave packet is changed by a 
given triplet of other wave packet modes \± — {fii, fX2, 1^3} 
as 



\4 



1)1 



fi 



^Vl n iJ,2 

Ry.u 



where 



Ru.fl 



E,, + E,, 



E„ 



En 



(9) 



(10) 



From now on, we assume that all the modes that belong 
to the packet (i.e., that are located between the two ex- 
ponential tails of the wave packet) have the same norm 
density equal to n. The perturbation approach breaks 



down and resonance sets in when «Jri < 



Substi 



tuting Eq. © for 

as 



one can rewrite the last inequality 



R v ,p, < fin. 



(11) 



This expression tells us that the resonance condition, for 
a given normal mode v, is fulfilled if there is at least one 
triplet of modes jl that satisfies inequality (jlll) . 

The probability for the onset of a resonance can there- 
fore be calculated with the following statistical numerical 
analysis @, [32J. For a given normal mode v, we define 
Ru,im, = xnha.fi{R Vj ft}. Collecting -R„ jA r for many modes 
and many values of the phase ip, we find the probability 
density distribution W(R vl r ). From this quantity we 
can calculate the probability V for a mode, with norm 
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FIG. 3: (color online) Comparison between the probability 
density function W(-R„ i( ij) of the quasiperiodic DNLS model 
and of the random DNLS model. For the quasiperiodic case, 
A = 2.5, while for the random case, we choose a disorder 
strength that gives a similar localization length. 



density n, to be resonant with at least one triplet of other 
modes at a given value of the interaction parameter fi. 
This is obtained by integrating W(R V n? ) from zero to 
fin 



V 



W{R) dR. 



(12) 



An example of probability density W(i?„ iA r ) for A = 2.5 
is shown in Fig. [3] (red line). For comparison we also 
show the same quantity for the random DNLS model 
(black line), as discussed in [f||3^|. Except for fine struc- 
tures, like small sharp peaks appearing in the quasiperi- 
odic case, the overall behaviour is qualitatively very simi- 
lar in the two cases. In particular, the probability density 
W tends to a finite constant value C when i?„ iA r — > 0. 
As a consequence, for small values of fin, a non-zero frac- 
tion of modes in the packet is resonant. The probability 
to be resonant is given by V ~ Cfin, thus we are in the 
weak chaos regime. For large values of fin, instead all 
the modes interact resonantly and V = 1; we are then in 
the strong chaos regime. 

Following the reasoning presented in this implies 
that also in the quasiperiodic case, as in disordered sys- 
tems, we may expect to find n%2 ~ t 1 / 3 in the weak chaos 
regime and 7712 ~ t 1 / 2 in the strong chaos regime. Note 
the strong chaos regime can only exist as a transient 
regime: as the wave packet spreads, its norm density 
n decreases, and eventually will reach a situation where 
fin < d. At this point, a crossover from strong to weak 
chaos is expected to occur during the time evolution [ll| . 

Let us finally stress that the "transition lines" that we 
have introduced by comparing the nonlinear frequency 
shift with the typical energy scales of the linear spec- 
trum do not define sharp phase transitions between dif- 
ferent spreading regimes. Instead, we may expect to see a 
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FIG. 4: (color online) Numerical results obtained by integrating the DNLS equations of motion ([2]). The time evolution of 
(log 10 7712} (left panel, top), 7 (right panel, top), (log 10 P) (left panel, bottom), and (£} (right panel, bottom) is shown versus 
log 10 t for different values of the nonlinear parameter j3 — 0.5, 1, 5, 10, 100. The initial wave packet in all simulations is a square 
distribution with L = 13 and the potential strength is A = 2.5. In the top right panel the two dashed lines correspond to 
theoretically predicted power laws 7 = 1/3 and 7 = 1/2. The width of the lines for the quantities (log 10 mz), (log 10 P) and 7 
represents the statistical error, which depends on time and on the number of realizations. In most cases the statistical error is 
smaller than the resolution of the figure. All quantities are dimensionless. 



relatively smooth crossover, such that the regimes of self- 
trapping, strong chaos and weak chaos should be clearly 
identified only far from the transition lines. 



V. TIME EVOLUTION 

We perform extensive numerical simulations solving 
Eqs. (J5J) and ([5]) for different sets of parameters {A, (3} 
and {A,£}, respectively. For each choice of parameters 
we average over N different realizations of the quasiperi- 
odic potential obtained by randomly changing the phase 
shift (p. For initial conditions, we use compact wave pack- 
ets that lie in the center of our computational box, taking 
care that during the time evolution the wave packet never 
reaches the box boundaries. The number of realizations 
considered varies between 100 and 500 and the number 
of lattice sites between 200 and 2000. To solve the equa- 
tions of motion, we use symplectic integration schemes 
of the SABA family 0, \$ that allow us to reach large 
integration times with good accuracy [36j . 



In order to quantify the type of subdiffusive behaviour, 
we calculate the exponent 7 by considering the logarithm 
of the second moment log 10 7712 for different realizations of 
the potential. We compute the average value (log 10 vn<i) 
and its statistical error, given by the standard deviation 
divided by the square root of the number of realizations 
N. Then the value of 7 at a given time t is calculated by 
applying a linear fitting procedure to the curve (log 10 7712) 
within a fixed time interval around log 10 t. By repeating 
this procedure at different t, we extract the behaviour of 
7 as a function of time and its relative statistical error. 



A. Results of the DNLS model 

Let us first show our results for the DNLS model. For 
the initial wave packet, we choose a square shaped distri- 
bution which equally populates L lattice sites with norm 
density rij = n = 1/L. In Fig. |||we present a represen- 
tative set of simulations for A = 2.5. We choose L = 13, 
which gives an initial localization volume larger than V. 
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The different panels show the time evolution of the sec- 
ond moment (log 10 m2), the spreading exponent 7, the 
participation ratio (log 10 P) , and the compactness index 
(£) . The width of the curves for (log 10 m-2) , (log 10 P) and 
7 corresponds to the statistical error. The values of the 
nonlinear frequency shift 5 induced by the initial wave 
packets used in these simulations are shown in Fig. [5] 
(empty downward triangles) in order to compare them 
to the relevant energy scales A and d. 

In all simulations we observe that nonlinearity causes 
the wave packet to spread. The spreading starts earlier 
when /3 is larger. We find that the spreading is always 
subdiffusive (7 < 1), confirming the result of previous 
works [13, • Subdiffusion is seen both in the second 
moment 7712 and in the participation ratio P, except for 
the largest value of /3 (yellow curves in Fig.HJ). In the lat- 
ter case, P saturates to a constant value after a transient 
time - a clear signature of self-trapping. This observa- 
tion of self-trapping only for j3 — 100 is consistent with 
the energy scale arguments schematically represented in 
Fig. [21 In the absence of self-trapping, the compactness 
index £ saturates to a constant value, indicating that the 
wave packet spreads but does not become more sparse. 
Conversely, in the presence of self-trapping the central 
part of the wave packet remains spatially trapped while 
its tails keep expanding, thus resulting in a wave packet 
that becomes more sparse during the evolution - nicely 
quantified by the compactness index which decreases to 
zero. We notice that the portion of packet that is ex- 
panding is characterized by a value of 7 larger than 1/2. 
After an initial increase, 7 reaches a maximum and then 
decreases to smaller values. In this regime, the evolution 
is rather complex. A similar behaviour was previously 
obtained also in random systems [ill . ll2T ] . The transient 
large values of 7 may be due to a nontrivial interacting 
mechanism that takes place between the expanding part 
and the self-trapped portion, resulting in faster spreading 
- an effect labeled as "overshooting" . 

For the lowest values of j3 the energy scale arguments 
suggest the occurrence of weak chaos. Indeed for (3 — 
0.5 and 1 the exponent saturates asymptotically around 
the theoretical value 7 = 1/3 (red and green curves in 
Fig. |4]), as expected. It is worth mentioning that this 
asymptotic exponent is the same as in random systems [3, 
111 ; meaning that the mechanism leading to destruction 
of exponential localization is rather universal. 

In difference to the random case, here during the time 
evolution, the value of 7 temporarily increases above 1/3, 
eventually reaching its asymptote only at longer times. 
This is an overshooting similar to the one that we have 
discussed above for the self-trapping regime, but occur- 
ring also for weaker nonlinearity. This effect is unique to 
the quasiperiodic system and is likely due to the presence 
of an infinite number of mini-bands and gaps in the linear 
spectrum of the Hamiltonian, which causes a temporary 
self-trapping of portions of the expanding wave packet in 
one or more energy gaps between mini-bands. This par- 
tial self-trapping is different from the self-trapping that 
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FIG. 5: (color online) Average logarithm of the second mo- 
ment of the expanding wave packet, (log 10 772,2} and spread- 
ing exponent, 7 for A = 2.2 (left plots) and A = 3.5 (right 
plots). For A = 2.2, the initial wave packet has width L — 31 
and we consider /3 = 0.18 (lower red curves), 1 (mid green 
curves) and 6.5 (upper blue curves). For A = 3.5, the ini- 
tial wave packet has width L — 5 and we consider /3 = 5.5 
(lower red curves), 15 (mid green curves), and 50 (upper blue 
curves). The width of the lines represents the statistical er- 
ror as in Fig. 0] Insets: average compactness index of the 
expanding wave packet (£) for the same sets of simulations. 



occurs when S > A, where all the packet modes are si- 
multaneously shifted out of resonance. For this reason 
partial self-trapping is not detectable as a saturation of 
the participation number P and can only be seen indi- 
rectly as an overshooting in the exponent 7. 

The two simulations for /3 = 5 and 10 lie in a range of 
energy were we expect to see strong chaos (blue and ma- 
genta curves in Fig. 2]) . As already said in the previous 
section, the strong chaos regime is transient: one should 
find a value of 7 around 1/2 for a few decades of time, 
eventually decreasing towards the asymptotic value 1/3. 
The two corresponding curves in Fig. 0] indeed exhibit a 
behaviour which qualitatively agrees with this expecta- 
tion. The value of 7 first rises up to 1/2, oscillates around 
this value and then starts to decrease as predicted. How- 
ever, especially for large (3, we also observe values of 7 
larger than 1/2. As in the weak chaos regime, this over- 
shooting again is evidence of partial self-trapping. Its 
mechanism is also transient and occurs in the same time 
intervals where strong chaos is expected. For this reason, 
while weak chaos is clearly observed in our simulations, 
strong chaos and partial self-trapping tend to overlap, 
thus producing a more complex evolution of the wave 
packet in quasiperiodic systems than in random systems. 

In Fig. [S]we show the results of simulations for A = 2.2 
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FIG. 6: (color online) Numerical results obtained by integrating the KG equations of motion ([5]). The time evolution of 
(log 10 7712} (left panel, top), 7 (right panel, top), (log 10 P) (left panel, bottom), and (£) (right panel, bottom) is shown versus 
log 10 t. The parameters are {X,S} = {2.5, 0.005} , {2.5, 0.01} , {2.5, 0.055} , {2.5, 0.075} , {2.5, 1.0}. We used an initial wave 
packet with width L = 13 for £ — 0.005, 0.01, 0.075, 1 and L = 11 for £ = 0.075. The width of the lines for the quantities 
(log 10 m2), (log 10 P) and 7 represents the statistical error as in Fig. [4] In the top right panel the two dashed lines correspond 
to theoretically predicted power laws 7 = 1/3 and 7 = 1/2. 



and A = 3.5; the corresponding values of nonlinear fre- 
quency shift are reported as triangles in Fig. [21 The 
values of L are L = 31 for A = 2.2 and L = 5 for 
A = 3.5, both larger than V. For {\,f3} = {2.2,0.18} 
and {A,/3} = {3.5,5.5} energy scale arguments predict 
weak chaos. We indeed find a spreading exponent which 
approaches asymptotically the value 1/3. For {A,/3} = 
{2.2,1}, {A,/3} = {2.2,6.5} and {A,/3} = {3.5,15} the 
predicted behaviour is either strong chaos or a regime in 
between strong and weak chaos. What we observe nu- 
merically is a growth of the spreading exponent 7 up to 
1/2 and even to larger values, followed by a decrease to- 
wards 1/3. In most cases, our simulations show a signifi- 
cant overshooting due to partial self-trapping. It is worth 
mentioning that this effect is larger for weaker disorder 
strength A, consistent with the fact the linear spectrum 
exhibits larger mini-gaps in this regime (see Fig. [1} . Fi- 
nally for {A,/3} = {3.5,50}, we observe self-trapping, as 
expected. 

In conclusion, from the analysis of the results of the 
DNLS model for different values of A we find that the en- 
ergy scale arguments and the model discussed in Section 
IV correctly explain the overall trend of the numerical 



simulations and the separation between different spread- 
ing regimes in the parameter space. 



B. Results of the KG model 

Due to the existence of a mapping between KG and 
DNLS, we expect to observe the same spreading regimes 
in the two models. This has been already proven in 
purely random systems where the two models reveal sim- 
ilar qualitative results in a wide range of parameters 
[1, 0, H, HH, EH- Despite this similarity, the study of 
the KG model remains interesting for at least two rea- 
sons. On one hand, it allows for testing the generality 
of the result in a case where there is just one conserved 
quantity. This is highly nontrivial - especially for self- 
trapping - for which rigorous results have been recently 
derived only in the case of Hamiltonians conserving both 
energy and norm Q. On the other hand, the KG model 
is advantageous from a numerical point of view. The fact 
that there is just one conserved quantity results in two 
orders of magnitude faster integration speed within the 
same integration error. 
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Similarly to what was done for the DNLS model, we 
initially set the compact wave packets to span a width 
L = 13 (unless otherwise stated) centered in the lat- 
tice, such that each site has equal energy £j = £ = 
%/L. This is implemented by setting initial momenta 
of p = ±V2£ with randomly assigned signs and zero co- 
ordinates. The values of initial energy densities £ are 
chosen to give expected spreading regimes of asymptotic 
weak chaos, intermediate strong chaos, and dynamical 
crossover from strong chaos to the slower weak chaos sub- 
diffusive spreading [9j. 

The results of the time simulations are shown in Fig. [51 
while the expected spreading regimes are given in Fig. [2] 
(full upward triangles) 37]. As one can see by com- 
paring Fig. [6J with Fig. [H the qualitative behaviour of 
the two models is rather similar. After initial transients, 
which increase with decreasing nonlinearity, all KG sim- 
ulations reveal subdiffusive growth of the second moment 
m.2 according to power law mi ~ t 7 with 7 < 1. If self- 
trapping is avoided, all simulations show a similar subdif- 
fusive behaviour for the participation numbers; moreover, 
the wave packets remain compact as they spread, since 
compactness indices at the largest computational times 
saturate around a constant (£) 3.5 ± 0.25. For the 
two smallest values of initial energy density £ — 0.05 
and £ — 0.01, the characteristics of the weak chaos 
regime are observed, namely, the exponent 7 saturates 
around 1/3 (red and green curves in Fig. \§§ after a tran- 
sient time. We stress that the only difference from the 
purely random systems is the overshooting phenomenon 
at transient times. This effect is an inherent property 
of quasiperiodic systems which inevitably manifests it- 
self in all spreading regimes, while in the disordered case 
it was shown to occur only in the regime of self-trapping 
00. 

For the two energy densities £ = 0.055 and 0.075 we 
suggest strong chaos, with characteristics similar to the 
DNLS case. The simulation with £ = 0.055 (blue curves 
in Fig. [6j indeed exhibits the typical behaviour of the 
strong chaos scenario: the characteristic exponent 7 in- 
creases up to predicted value 1/2 and remains so for 
about two time decades, followed by a crossover with 
7 decreasing to the weak chaos dynamics. There is also 
another possibility for larger £ = 0.075, when interme- 
diate strong chaos is masked due to partial self-trapping 
(magenta curves in Fig. [6]). Thus, 7 shows values larger 
then 1/2 but still with subsequent decay to slower subd- 
iffusion. Here, we would like to strongly emphasize that 
none of the simulations exhibit pronounced deviations 
from strong or weak chaos regimes of spreading, i.e. long- 
lasting overshooting with 7 > 1/2, or significant slowing 
down to values 7 < 1/3. 

Finally, for £ = 1.0 the dynamics enter the self- 
trapping regime, as our theory predicts. There a major 
part of the initial excitation stays localized, while the 
remainder spreads (yellow curves in Fig. [S]). The par- 
ticipation numbers, therefore, do not grow significantly 
and (log 10 P) starts to level off at large time (Fig. [51 



left panel, bottom, yellow curve). In contrast, the small 
spreading portion yields a continuous increase of the sec- 
ond moment m 2 (Fig. [5J left panel, top, yellow curve), 
which initially is characterized by large values of 7 > 1/2 
(howbeit, for larger time 7 decreases). Consequently, the 
compactness index (£) (Fig. [Bl right panel, bottom, yel- 
low curve) drops down to small values indicating deep 
self-trapping regime. Note that a similar behaviour has 
been observed before in purely random systems [IJ (HI ■ 
Unusually large values of mi can be explained by lo- 
cal trapping-detrapping processes in the evolving wave 
packet. The corresponding dynamics is in strong non- 
equilibrium - its theoretical description has yet to be 
developed. 



C. Role of the shape of the initial wave packet 

In this subsection we show that the results discussed 
so far do not depend on the shape of the initial wave 
packet. Besides its theoretical interest, this issue is also 
relevant from the point of view of experiments, where it 
is not always possible to design the wave packets at will. 

In the previous sections, we have used a square dis- 
tribution as the initial wavepacket. Now, inspired by 
the experiments with ultracold atoms, we consider initial 
wave packets with the shape of a Gaussian distribution or 
a Thomas-Fermi (TF) distribution. The Gaussian wave 
packet centered around the site J = has the form 

^■(0) = Cie-£, (13) 

where a is a parameter controlling the width of the packet 
while C\ is a constant factor that can be determined 
by using the normalization condition Ylj |^-| = 1. A 
Thomas-Fermi wave packet is instead defined by 

^(0)=C 2] jl^ (14) 

in the region where \j\ < R and ipj = otherwise. The 
parameter R is the Thomas-Fermi radius characterizing 
the width of the distribution, while the constant Ci is a 
normalization factor. These two distributions are of in- 
terest when considering ultracold bosons initially released 
from an harmonic trap in the Gross-Pitaevskii regime 

In Fig. [7] we show the time evolution in the DNLS 
model of the second moment of the expanding wave 
packet, (log 10 77i2) (top row) and of the spreading ex- 
ponent, 7 (bottom row), using initially a Gaussian (left 
column) and a TF (right column) wave packet distribu- 
tion. In the insets we also show the compactness index 
(£), in order to identify the self-trapping regime. We 
choose the width of the initial distributions (er and R) so 
that the nonlinear frequency shift is similar to the one 
already used for the simulations in Fig. 0] 39]. In par- 
ticular we use cr = 5 and R = 7.5, yielding a nonlinear 
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FIG. 7: (color online) Average logarithm of the second mo- 
ment of the expanding wave packet, (log 10 7712} and spreading 
exponent, 7 as a function of time for different nonlinearities 
/3 = 0.5, 1, 5, 10, 100. The disorder strength is A = 2.5 in all 
simulations. As an initial condition, we have used a Gaussian 
wave packet with cr = 5 (left plots) and a TF distribution 
with R — 7.50 (right plots). The width of the lines represents 
the statistical error as in Fig. [4] Insets: average compactness 
index of the expanding wave packet (£} for the same sets of 
simulations. 



frequency shift S « /3/13. The values of (3 used in Fig. [7] 
are the same as those previously considered. 

From the comparison between the results of Fig. [7] and 
Fig. [H we can conclude that the shape of the initial wave 
packet does not affect the overall behaviour of the time 
evolution, nor its interpretation in terms of regimes of 
weak and strong chaos, self-trapping, and overshooting. 
This suggests the results that we have obtained are rather 
general and that the nonlinear frequency shift S is the 
only key parameter controlling the dynamics of the wave 
packet. 



D. Application to cold atoms 

The DNLS model can be used to simulate the dynam- 
ics of bosons in optical lattices at zero temperature [24| 
and in the tight-binding regime, where the DNLS equa- 
tion corresponds to a discretized version of the Gross- 
Pitaevskii (GP) equation for the dynamics of a Bosc- 
Einstein condensate in the single-band approximation. 
The validity of this mean-field theory is not ensured for 
those dynamical regimes where GP predicts chaos [40| . 
which can be viewed as a signature of a large depletion 
of the condensate. For this reason, in the presence of dis- 
order the theory fails to predict the long time evolution 



of observables directly related to small scale fluctuations 
and long-range coherence. However, for coarse-grained 
observables, like the width of the wave packet in real and 
momentum space, or the participation number, the pre- 
dictions of the theory remain very good even in regimes 
where the depletion is expected to be large, long after the 
random fluctuations prevent the prediction of fine scale 
structures. This has been recently shown in Ref. 41| by 



comparing the predictions of the GP equation with one 
beyond mean-field theory in numerical simulations within 
timescales of the order of typical experiments with cold 
atoms and long enough to observe the effects of deple- 
tion and chaotic dynamics. Indeed our analysis is essen- 
tially based on coarse-grained observables. In addition, 
for each set of parameters we also average over many 
realizations and this extends the validity of the present 
approach even for longer times, as any residual depen- 
dence on small scale fluctuations is further suppressed 
by the averaging procedure. 

When applied to bosons expanding in bichromatic op- 
tical lattices, our results provide a consistent interpreta- 
tion of the experimental data of Ref. [22J , where a Bose- 
Einstein condensate, initially confined in an harmonic 
trap, is let free to expand in a bichromatic potential. In 
our dimensionless units, the expansion lasts for times of 
the order of 10 4 (see [20| for more details) and the width 
of the atomic cloud increases up to 50 — 100 lattice sites. 
In this experiment a subdiffusivc spreading is observed 
with exponents 7 significantly larger than 1 /3 already for 
weak nonlinearities and even larger than 1/2 for larger 
nonlinearities [421 ] . Our work suggests that such large 
values of 7 can be explained in terms of a transient over- 
shooting caused by partial self-trapping in mini-bands. 



VI. SUMMARY AND CONCLUSIONS 

In this work we have considered the problem of the 
interplay between localization and interaction in one- 
dimensional quasiperiodic systems. We have investigated 
the expansion of initially localized wave packets in two 
different quasiperiodic models, DNLS and KG. We have 
confirmed that interaction destroys localization, giving 
rise to a subdiffusivc growth of the second moment of 
the wave packet (rri2 ~ t 1 with 7 < 1). 

We have interpreted the spreading process in terms 
of resonances in the mode-mode coupling. In particu- 
lar, we have identified the different spreading regimes of 
self-trapping, strong chaos and weak chaos by comparing 
the frequency shift induced by the nonlinearity with the 
energy scales extracted from the spectrum of the under- 
lying linear system. For weak and strong chaos regimes 
we have also predicted the expected spreading exponents 
7=1/3 and 7=1/2 respectively. 

We have performed numerical simulations, which last 
for much longer times than existing simulations, and we 
have averaged our results over many realizations. This 
gave us the possibility to accurately calculate the spread- 
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ing exponent 7 and observe the weak chaos regime. A 
key difference with respect to random systems @, Q is 
the occurrence of transient overshooting regimes that we 
interpret as due to the peculiar structure of the linear 
spectrum of the quasiperiodic system, which is separated 
into mini-bands. These mini-bands are responsible for 
mechanisms of partial self-trapping. Signatures of strong 
chaos have also been observed, but the temporal over- 
lap of strong chaos and partial self-trapping makes the 
analysis of the spreading more complex than for random 
systems. We have also verified that our main results do 
not depend on the details of the shape of the initial wave 
packet. This suggests that the nonlinear frequency shift, 
(5, is the only parameter that controls the dynamics. 
Finally, our results provide a consistent interpretation 



of the subdiffusive spreading observed in experiments 
with ultracold atoms propagating in bichromatic optical 
lattices [22l |. 



Acknowledgments 

M.L. thanks the Max Planck Institute for the Physics 
of Complex Systems, Dresden, for their hospitality. We 
are indebted to G. Modugno and E. Lucioni for fruit- 
ful discussions. This work has been supported by ERC 
through the QGBE grant and by the UPV/EHU under 
program UFI 11/55. 



[1] P. W. Anderson, Phys. Rev. 109, 1492 (1958). 
[2] B. Kramer, A. MacKinnon, Rep. Prog. Phys. 56, 1469 
(1993). 

[3] G. Kopidakis, S. Komineas, S. Flach, and S. Aubry, Phys. 

Rev. Lett. 100, 084103 (2008). 
[4] A. S. Pikovsky and D. L. Shepelyansky, Phys. Rev. Lett. 

100, 094101 (2008). 
[5] I. Garcfa-Mata and D. L. Shepelyansky, Phys. Rev. E 79, 

026205 (2009). 

[6] S. Flach, D. O. Krimer, and Ch. Skokos, Phys. Rev. Lett. 

102, 024101 (2009). 
[7] Ch. Skokos, D. O. Krimer, S. Komineas, and S. Flach, 

Phys. Rev. E 79, 056211 (2009). 
[8] H. Veksler, Y. Krivolapov, and S. Fishman, Phys. Rev. 

E 80, 037201 (2009) 
[9] S. Flach, Chem. Phys. 375, 548 (2010). 
[10] C. Skokos and S. Flach, Phys. Rev. E 82, 016208 (2010). 
[11] T. V. Laptyeva, J. D. Bodyfelt, D. O. Krimer, C. Skokos, 

and S. Flach, Europhys. Lett. 91, 30001 (2010). 
[12] J. D. Bodyfelt, T. V. Laptyeva, C. Skokos, D. O. Krimer, 

and S. Flach, Phys. Rev. E 84, 016205 (2011). 
[13] J. D. Bodyfelt, T. Kottos, and B. Shapiro, Phys. Rev. 

Lett. 104, 164102 (2010). 
[14] H. Hiramoto and M. Kohmoto, Int. J. Mod. Phys. B 6, 
281 (1992). 

[15] S. Aubry and G. Andre, Ann. Israel Phys. Soc. 3, 133 
(1980). 

[16] Y. Lahini, R. Pugatch, F. Pozzi, M. Sorel, R. Morandotti, 
N. Davidson, and Y. Silberberg, Phys. Rev. Lett. 103, 
013901 (2009). 

[17] G. Roati, C. D'Errico, L. Fallani, M. Fattori, C. Fort, M. 

Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, 

Nature 453, 895 (2008). 
[18] G.S. Ng and T. Kottos, Phys. Rev. B 75, 205120 (2007). 
[19] M. Johansson, M. Hornquist, and R. Riklund, Phys. Rev. 

B 52, 231 (1995). 
[20] M. Larcher, F. Dalfovo, and M. Modugno, Phys. Rev. A 

80, 053606 (2009). 
[21] Z. Zhang, P. Tong, J. Gong, and B. Li, Phys. Rev. E 83, 

056205 (20011). 
[22] E. Lucioni, B. Deissler, L. Tanzi, G. Roati, M. Zaccanti, 

M. Modugno, M. Larcher, F. Dalfovo, M. Inguscio, and 

G. Modugno, Phys. Rev. Lett. 106, 230403 (2011). 



[23] We note that our choice of a is equivalent to the more 
commonly used value of the inverse of the golden mean 

(VE-i)/2. 

[24] A. Trombettoni and A.Smerzi, Phys. Rev. Lett. 86, 2353 
(2001). 

[25] M. Modugno, New J. Phys. 11, 033023 (2009). 

[26] D. N. Christodoulides, F. Lederer, and Y. Silberberg, 
Nature 424, 817 (2003). 

[27] A.A. Ovchinnikov, N.S. Erikhman and K.A. Pronin, 
Vibrational-Rotational Excitations in Nonlinear Molecu- 
lar Systems, Kluwer Academic/Plenum Publishers, New 
York (2001). 

[28] Yu. S. Kivshar, M. Peyrard, Phys. Rev. A 46, 3198 
(1992); Yu. S. Kivshar, Phys. Lett. A 173, 72 (1993); 
M. Johansson, Physica D 216, 62 (2006). 

[29] H. Hiramoto and S. Abe, J. Phys. Soc. Japan 57, 1365 
(1988). 

[30] S. Y. Jitomirskaya, Ann. Math 150, 1159 (1999). 

[31] M. Larcher, M. Modugno, and F. Dalfovo, Phys. Rev. A 

83, 013624 (2011). 
[32] D. O. Krimer and S. Flach, Phys. Rev. E 82, 046221 

(2010). 

[33] An intuitive understanding of this band structure can be 
given, following an heuristic argument. The wavelength 
of the potential Vj is l/|a| = 2.618 . . . An effective wave- 
length equal to an integer number q would correspond to 
a separation in exactly q bands. Our value of l/|a| lies 
between two and three, so that the band structure has 
neither two nor three bands, but three main bands with 
an internal structure of sub-bands. 

[34] The self-trapping regime for the DNLS case can be un- 
derstood also in terms of energy and particle conserva- 
tion [3, 0, E3] • Due to particle conservation, the particle 
density decreases during the spreading, which implies a 
decreasing of the mean-field interaction energy. If S > A, 
the nonlinearity is so dominant that the linear part of the 
Hamiltonian is unable to compensate the loss of mean- 
field energy; therefore a part of the wave packet does 
not spread and stays localized. The presence of gaps in 
the linear spectrum is then a necessary condition for the 
occurrence of self-trapping. 

[35] J. Laskar and P. Robutel, Celest. Mech. Dyn. Astron. 80, 
39 (2001). 



12 



[36] The numerical accuracy of our calculation is controlled by 
checking the conservation of the energy H and the norm 
S of the expanding wave packet. The error is always kept 
smaller than 10 -2 ' 5 . For the integration we used time 
steps between 0.1 and 0.05. 

[37] In order to compare the nonlinear frequency shift 5 ~ £ of 
the KG model with the energy scales of the DNLS model 
A and d, we use the approximate mapping j3n ~ 6A£ at 
t = 0. Therefore, the quantity that is plotted in Fig. [2] 
for the KG model is 6A£ . 

[38] F. Dalfovo, S. Giorgini, L.P. Pitaesvkii, and S. Stringari, 
Rev. Mod. Phys. 71, 463 (1999). 

[39] For Gaussian and TF distributions we identify the non- 
linear frequency shift with the quantity /3 . j 4 , which 
is also identical (up to a prefactor) with the mean-field 
interaction energy in Eq. (1). 



[40] Y. Castin, and R. Dum, Phys. Rev. Lett. 79, 3553 
(1997); S.A. Gardiner, D. Jaksch, R. Dum, J.I. Cirac, 
and P. Zoller, Phys. Rev. A 62, 023612 (2000); C. Weiss, 
and N. Teichmann, Phys. Rev. Lett. 100, 140408 (2008); 
J. Reslen, C.E. Creffield, and T.S. Monteiro, Phys. Rev. 
A 77, 043621 (2008); I. Bfezinova, L.A. Collins, K. Lud- 
wig, B.I. Schneider, and J. Burgdorfer, Phys. Rev. A 83, 
043611 (2011); M. Heimsoth, C.E. Creffield, L.D. Carr, 
and F. Sols, New. J. Phys. 14, 075023 (2012). 

[41] I. Bfezinova, A. Lode, A. Streltsov, O. Alon, L. Ceder- 
baum, and J. Burgdorfer, Phys. Rev. A 86, 013630 
(2012). 

[42] We notice that the spreading exponent a in Ref. [22] 
differs from our exponent 7 by a factor two, a = 7/2, 
due to a different notation. 



